Faddeev_Leverrier Subroutine

public subroutine Faddeev_Leverrier(A, c, Ainv, success, check)

Arguments

Type IntentOptional Attributes Name
real(kind=dp), intent(in), DIMENSION(:, :) :: A
real(kind=dp), intent(out), DIMENSION(:) :: c
real(kind=dp), intent(out), optional, DIMENSION(SIZE(A,1), SIZE(A,1)) :: Ainv
logical, intent(out), optional :: success
logical, intent(in), optional :: check

Calls

proc~~faddeev_leverrier~~CallsGraph proc~faddeev_leverrier Faddeev_Leverrier proc~identity_n Identity_n proc~faddeev_leverrier->proc~identity_n proc~trace Trace proc~faddeev_leverrier->proc~trace

Variables

Type Visibility Attributes Name Initial
integer, public, parameter :: dp = KIND(1.0d0)
real(kind=dp), public, DIMENSION(SIZE(A,1), SIZE(A,1)) :: Bk
real(kind=dp), public, DIMENSION(SIZE(A,1), SIZE(A,1)) :: I
real(kind=dp), public, DIMENSION(SIZE(A,1), SIZE(A,1)) :: B_Nm1
real(kind=dp), public, DIMENSION(SIZE(A,1), SIZE(A,1)) :: AB
logical, public :: do_check = .TRUE.
integer, public :: N
integer, public :: k

Source Code

    SUBROUTINE Faddeev_Leverrier(A, c, Ainv, success, check)
        INTEGER, PARAMETER :: dp = KIND(1.0d0)
        REAL(dp), DIMENSION(:, :), INTENT(IN)  :: A
        LOGICAL, OPTIONAL, INTENT(IN) :: check
        REAL(dp), DIMENSION(:),     INTENT(OUT) :: c
        REAL(dp), DIMENSION(SIZE(A,1), SIZE(A,1)), OPTIONAL, INTENT(OUT) :: Ainv
        LOGICAL, OPTIONAL, INTENT(OUT) :: success
        REAL(dp), DIMENSION(SIZE(A,1), SIZE(A,1)) :: Bk, I, B_Nm1, AB
        LOGICAL :: do_check = .TRUE.
        INTEGER :: N, k

        N = SIZE(A,1)

        IF (PRESENT(check)) do_check = check

        IF (do_check) THEN
            PRINT*, "Checking if the matrix A is square and size of c is correct"
            IF (SIZE(A,2) /= N .OR. SIZE(c) < N+1) THEN
                PRINT *, "Error : Matrix A must be square and size of c must be at least N+1"
                STOP
            END IF
        END IF

        ! Initialization
        I = Identity_n(N)
        c = 0.0_dp
        c(1) = 1.0_dp
        c(2) = -Trace(A)
        Bk = A + c(2)*I

        DO k = 2, N
            AB = MATMUL(A, Bk)
            c(k+1) = -Trace(AB) / REAL(k, dp)
            Bk = AB + c(k+1)*I
            IF (k == N-1 .AND. PRESENT(Ainv)) B_Nm1 = -Bk
        END DO

        IF (PRESENT(Ainv) .AND. PRESENT(success)) THEN
            IF (ABS(c(N+1)) < 1.0e-12_dp) THEN
                success = .FALSE.
                Ainv = 0.0_dp
            ELSE
                success = .TRUE.
                Ainv = B_Nm1 / c(N+1)
            END IF
        ELSE IF (PRESENT(Ainv)) THEN
            IF (ABS(c(N+1)) < 1.0e-12_dp) THEN
                Ainv = 0.0_dp
            ELSE
                Ainv = B_Nm1 / c(N+1)
            END IF
        END IF

    END SUBROUTINE Faddeev_Leverrier